CalMAPPER

Description

This script merges CalMAPPER activity data with treatment polygons. This is a necessary step before analyzing prescribed fire data in CalMAPPER.

Define path to Calmapper and to the gdb file

calmapper_dir <- fs::dir_ls(ref_path, recurse = T, glob = '*CalMAPPER', type = 'directory')
calmapper_gdb <- fs::dir_ls(calmapper_dir, recurse = T, glob = '*.gdb')
if(length(calmapper_gdb) > 1){
  stop( "There's more than one CalMAPPER .gdb file. script assumes only 1.")
}

st_layers(calmapper_gdb)

## Driver: OpenFileGDB 
## Available layers:
##                 layer_name     geometry_type features fields
## 1 CMDash_ProjectTreatments     Multi Polygon     1560     21
## 2     CMDash_TreatmentPols     Multi Polygon     2784     23
## 3    CMDash_TreatmentLines Multi Line String       23     22
## 4     CMDash_TreatmentPnts       Multi Point      136     20
## 5        CMDash_Activities                NA    20439     38
## 6          CMDash_Metadata                NA        1      7
## 7       CalMapper_fire_act     Multi Polygon     1861     28
##                   crs_name
## 1 WGS 84 / Pseudo-Mercator
## 2 WGS 84 / Pseudo-Mercator
## 3 WGS 84 / Pseudo-Mercator
## 4 WGS 84 / Pseudo-Mercator
## 5                     <NA>
## 6                     <NA>
## 7 WGS 84 / Pseudo-Mercator

act <- st_read(calmapper_gdb, layer = 'CMDash_Activities')

## Reading layer `CMDash_Activities' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb' 
##   using driver `OpenFileGDB'

## Warning: no simple feature geometries present: returning a data.frame or tbl_df

trt <- st_read(calmapper_gdb, layer = 'CMDash_TreatmentPols')

## Reading layer `CMDash_TreatmentPols' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CalMAPPER\CALFIRE_FuelReductionProjects_2023\CALFIRE_FuelReductionProjects.gdb' 
##   using driver `OpenFileGDB'

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts. The
## processing may be really slow.  You can skip the processing by setting
## METHOD=SKIP, or only make it analyze counter-clock wise parts by setting
## METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock
## wise defined

## Simple feature collection with 2784 features and 23 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13842530 ymin: 3842624 xmax: -12941530 ymax: 5161286
## Projected CRS: WGS 84 / Pseudo-Mercator

Clean

Remove non-treatment “activities”

act <- act %>% 
  filter(!ACTIVITY_DESCRIPTION %in% c("GIS Validation", "Education Outreach", "Project Administration", "Planning Meeting", "Public Contacts", "Water Site Development"))

Remove activities that aren’t complete or aren’t measured in acres

act <- act %>% 
  filter(ACTIVITY_STATUS == "Complete", UNIT_OF_MEASURE == "Acres")

Remove unnecessary columns

act <- act %>% 
  select(-UNIT_OF_MEASURE, -PROJECT_STATUS, -QUANTITY, -ACTIVITY_STATUS, -GROUND_DISTURBING, -ThisFiscal, -LastFiscal, -PrevFiscals, -OBJECTID_1, -CFIP_FR, -CONTRACT)

Merge Rx fire Activity data with treatment polygons

act_sf <- right_join(trt %>% select(PROJECT_ID, TREATMENT_ID), 
           act)

## Joining with `by = join_by(PROJECT_ID, TREATMENT_ID)`

Filter to fire

Find relevant treatment names

act_sf %>% 
  st_drop_geometry() %>% 
  group_by(ACTIVITY_DESCRIPTION) %>% 
  count() %>% 
  print(n=50)

## # A tibble: 29 × 2
## # Groups:   ACTIVITY_DESCRIPTION [29]
##    ACTIVITY_DESCRIPTION                      n
##    <chr>                                 <int>
##  1 Broadcast Burn                          545
##  2 Chaining                                 57
##  3 Chipping                               1701
##  4 Commercial Thinning (Cable Yarding)      18
##  5 Commercial Thinning (Tractor Yarding)    24
##  6 Crushing                                 35
##  7 Erosion Control                           5
##  8 Follow up - Herbicide                    56
##  9 Follow up - Other                         6
## 10 Follow up - Slash disposal              170
## 11 Fuel Break (Shaded)                      25
## 12 Grazing                                  23
## 13 Herbicide (Post-Treatment)               26
## 14 Herbicide (Pre-Treatment)                 8
## 15 Land Conservation                         4
## 16 Lop and Scatter                         521
## 17 Mastication                             766
## 18 Pile Burning                           1316
## 19 Piling (Manual)                        1919
## 20 Piling (Mechanical)                     251
## 21 Pruning                                 305
## 22 Rangeland Mowing                         55
## 23 Release - Herbicide                      15
## 24 Release - Mechanical                     93
## 25 Release - Other                          10
## 26 Site Preparation (CFIP)                  68
## 27 Thinning                                285
## 28 Thinning (Manual)                      2903
## 29 Thinning (Mechanical)                   200

Filter by treatment

act_sf_fire <- 
  act_sf %>% 
  filter(ACTIVITY_DESCRIPTION %in% c("Broadcast Burn", "Cultural Burning", "Pile Burning"))

Clean up dates

Remove hour information

test <- act_sf_fire %>% 
  mutate(ACTIVITY_START_test = as.Date(ACTIVITY_START)) %>% 
  mutate(ACTIVITY_END_test = as.Date(ACTIVITY_END)) 

Check

This should be 2/11 - 2/11. The values that R automatically reads in from ArcGIS are one day off.

test %>% 
  filter(AQ_ID == 69461) %>% 
  select(ACTIVITY_START, ACTIVITY_START_test, ACTIVITY_END, ACTIVITY_END_test)

## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13709070 ymin: 4748703 xmax: -13705400 ymax: 4751503
## Projected CRS: WGS 84 / Pseudo-Mercator
##        ACTIVITY_START ACTIVITY_START_test        ACTIVITY_END ACTIVITY_END_test
## 1 2020-02-10 16:00:00          2020-02-11 2020-02-10 16:00:00        2020-02-11
##                            SHAPE
## 1 MULTIPOLYGON (((-13706330 4...

Apply to everything

act_sf_fire <- act_sf_fire %>% 
  mutate(ACTIVITY_START = as.Date(ACTIVITY_START)) %>% 
  mutate(ACTIVITY_END = as.Date(ACTIVITY_END)) 

Create DURATION and YEAR columns

DURATION

act_sf_fire <- act_sf_fire %>% 
  mutate(DURATION = ACTIVITY_END - ACTIVITY_START+1)

Add YEAR column

act_sf_fire <- act_sf_fire %>% 
  mutate(YEAR = year(ACTIVITY_END))

Check

act_sf_fire %>% 
  filter(AQ_ID == 69461) %>% 
  select(ACTIVITY_START, ACTIVITY_END, DURATION, YEAR) 

## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13709070 ymin: 4748703 xmax: -13705400 ymax: 4751503
## Projected CRS: WGS 84 / Pseudo-Mercator
##   ACTIVITY_START ACTIVITY_END DURATION YEAR                          SHAPE
## 1     2020-02-11   2020-02-11   1 days 2020 MULTIPOLYGON (((-13706330 4...

Save data from all years, all places

cm <- act_sf_fire
save(cm, file = here("Rdata/CalMapper_activities_fire.Rdata"))
write.csv(cm, file = here("excel/CalMapper_activities_fire.csv"), row.names = F)

write.csv(act_sf_fire %>% st_drop_geometry(), file = "~/Reference data/CalMapper/activities_fire.csv", row.names = F)
write.csv(trt %>% st_drop_geometry(), file = "~/Reference data/CalMapper/treatments_fire.csv", row.names = F)

Are there Rx fires recorded in Treatments data but not in Activities data or vice versa?

Create subsets of data for completed broadcast burns

Activities

act_broadcast <- act_sf_fire %>% 
  filter(ACTIVITY_DESCRIPTION == "Broadcast Burn")

Treatments

trt %>% 
  st_drop_geometry() %>% 
  group_by(TREATMENT_OBJECTIVE) %>% 
  count()

## # A tibble: 5 × 2
## # Groups:   TREATMENT_OBJECTIVE [5]
##   TREATMENT_OBJECTIVE        n
##   <chr>                  <int>
## 1 Broadcast Burn           559
## 2 Forestland Stewardship   296
## 3 Fuel Break               215
## 4 Fuel Reduction          1542
## 5 Right of Way Clearance   172

trt_broadcast_complete <- trt %>% 
  filter(TREATMENT_OBJECTIVE == "Broadcast Burn") %>% 
  filter(ACTIVITY_STATUS == "Complete")

Check if there are treatment IDs in act_broadcast that aren’t in trt

check <- act_broadcast %>% 
  filter(!TREATMENT_ID %in% trt$TREATMENT_ID) %>% 
  nrow()

check2 <- act_broadcast %>% 
  filter(!TREATMENT_NAME %in% trt$TREATMENT_NAME) %>% 
  nrow()

if(check==0 & check2 == 0){
  paste("There are no TREATMENT_ID values or TREATMENT_NAME values in activities that aren't in trt, at least for broadcast burns")
} else{
  print("Alert! There are TREATMENT_ID or TREATMENT_NAME values in activities data that aren't in trt")
}

## [1] "There are no TREATMENT_ID values or TREATMENT_NAME values in activities that aren't in trt, at least for broadcast burns"

Check if there are treatment IDs in completed broadcast treatments that aren’t in activities

not_in_act <- trt_broadcast_complete %>% 
  filter(!TREATMENT_NAME %in% act_broadcast$TREATMENT_NAME) %>% 
  nrow()

print(paste("Alert! There are", not_in_act, "records in trt_broadcast_complete that aren't in activities. These are potential missing data from the final output!"))

## [1] "Alert! There are 7 records in trt_broadcast_complete that aren't in activities. These are potential missing data from the final output!"

Export those to put in .ppt

missing_from_activities <- trt_broadcast_complete %>% 
  filter(!TREATMENT_ID %in% act_broadcast$TREATMENT_ID) %>% 
  select(CALMAPPER_ID, PROJECT_NAME, TREATMENT_NAME, TREATMENT_OBJECTIVE, PROJECT_STATUS, ACTIVITY_STATUS, PROJECT_START_DATE, PROJECT_END_DATE, TREATMENTAREA_ACRES) 
missing_from_activities %>% 
  st_drop_geometry() %>% 
  summarize(missing_acres = sum(TREATMENTAREA_ACRES))

##   missing_acres
## 1        653.99

write_excel_csv(missing_from_activities %>% st_drop_geometry(), "calmapper_trt_not_act.xls")

Visualize broadcast treatments with no overlapping activities

map <- mapview(list(trt_broadcast_complete, act_broadcast, missing_from_activities), col.regions=list("red","blue", "green"),col=list("red","blue", "green"))
map

It’s safe to use the Activities information and not the Treatment information, aside from those 7 missing records.

Write activities data to gdb

st_write(act_sf_fire, dsn = calmapper_gdb, layer = "CalMapper_fire_act", delete_layer = T)

## Warning in clean_columns(as.data.frame(obj), factorsAsCharacter): Dropping
## column(s) DURATION of class(es) difftime

## Deleting layer `CalMapper_fire_act' using driver `OpenFileGDB'
## Writing layer `CalMapper_fire_act' to data source 
##   `C:/Users/ctubbesi/OneDrive - California Air Resources Board/Documents/Reference data/CalMAPPER/CALFIRE_FuelReductionProjects_2023/CALFIRE_FuelReductionProjects.gdb' using driver `OpenFileGDB'
## Writing 1861 features with 28 fields and geometry type Multi Polygon.

Find GIS acreage of “Upper Little Stony Post Ranch”

area_m2 <- act_sf_fire %>% 
  filter(AQ_ID == 267909) %>% 
  st_area()
area_ac <- area_m2/4046.86
area_ac

##  [m^2]

Save 2020 to use for Venn Diagram

act_2020 <- act_sf_fire %>% 
  filter(YEAR == 2020)

Mask out federal lands

Load NF data

load(file = here("Rdata/NF.Rdata"))

Clip

act_sf_fire_noNF <- act_sf_fire %>% 
  st_difference(NF)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

Save data from all years, outside NF

cm <- act_sf_fire_noNF
save(cm, file = here("Rdata/CalMapper_activities_fire_noNF.Rdata"))

act_2020 <- st_difference(act_2020, NF)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

Save

save(act_2020, file = here("Rdata/CalMapper_activities_fire_2020_noNF.Rdata"))